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We consider aspects of the population dynamics, inside a bound domain, of diffusing agents 
carrying an attribute which is stochastically destroyed upon contact with the boundary. The nor- 
mal mode analysis of the relevant Helmholtz equation under the partially absorbing, but uniform, 
boundary condition provides a starting framework in understanding detailed evolution dynamics of 
the attribute in the time domain. In particular, the boundary-localized depletion has been widely 
employed in practical applications that depend on geometry of various porous media such as rocks, 
cement, bones, and cheese. While direct relationship between the pore geometry and the diffusion- 
relaxation spectrum forms the basis for such applications and has been extensively studied, relatively 
less attention has been paid to the spatial variation of the boundary condition. In this work, we 
focus on the way the pore geometry and the inhomogeneous depletion strength of the boundary be- 
come intertwined and thus obscure the direct relationship between the spectrum and the geometry. 
It is often impossible to gauge experimentally the degree to which such interference occur. We fill 
this gap by perturbatively incorporating classes of spatially-varying boundary conditions and derive 
their consequences that are observable through numerical simulations or controlled experiments on 
glass bead packs and artificially fabricated porous media. We identify features of the spectrum that 
are most sensitive to the inhomogeneity and apply the method to the spherical pore with a sim- 
ple hemi-spherical binary distribution of the depletion strength and obtain bounds for the induced 
change in the slowest relaxation mode. 

PACS numbers: 89.90.+n,76.60.-k ,81.05.Rm,91.60.-x 



I. INTRODUCTION 

We consider the evolution of a physical attribute car- 
ried by a population of random-walkers inside a medium 
bound by a wall of general shape. When a walker hits 
the boundary, the encounter depletes its attribute with a 
certain probability, p g [0, 1]. Our main concern here is 
on allowing this probability to have general spatial vari- 
ation and to investigate its consequences on the spatio- 
temporal evolution of the local attribute density ^'(r,t) 
and its net sum, M.{t) = J (ir^'(r,t). Without the spa- 
tial variations of p and the local difFusivity, D, the prob- 
lem reduces to the classic Helmholtz equation with a 
uniform Robin's boundary condition, bookended by the 
Dirchlet- (p — > oo) on one limit and by the Neumann- 
condition (p — !■ 0) on the other. The spectral analysis 
of its eigenmodes has been discussed as a probe of ge- 
ometrical properties of the boundary [TJ [51 [3] and found 
application for a variety of systems such as the electrode 
impedance [4|, acoustics 0, NMR relaxometry[6j |7], nu- 
clear level statistics [H], quantum chaos [S], and migration 
of cultural or genetic trait [101 El HI- Population evo- 
lution of the Web crawler-programs [13] deployed over a 
large network in the presence of unstable nodes may be 
an example where the diffusion may not necessarily be 
bound to the physical space. 

To be concrete, we are directly motivated by issues en- 
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countered in the interpretation of the magnetic resonance 
(MR) probe of fluid in conventional porous media[T4l [15] . 
suspended particulate aggregates and colloids [TB I 1^71 [T5] . 
The diffusion-relaxation dynamics of polarized proton 
spins carried by diffusing molecules has been widely ex- 
ploited, sometimes without full justification, to char- 
acterize the pore geometry and fluid viscosity in the 
soilTS', cements [H], oil pigment of old paintings [22], 
biological tissues [35], flber bundles [31]. plant cells [25 or 
a piece of cheese [26]. In its geophysical or oil-fleld ap- 
plication, the interface-enhanced relaxation is used as a 
probe for the pore geometry of rocks ^l5. ,27. .28, , composi- 
tion of pore flUing fluids [35] . and even wetting conditions. 
This extraordinary utility derives from the basic observa- 
tion that the interface-enhanced relaxation rate (widely 
called T2 distribution [TS]) is directly proportional to the 
surface-to- volume ratio of the pore enclosure when cer- 
tain conditions are met. (See Eqs T]|2 below) 

Three basic conditions are required to be met implic- 
itly in such a mapping between the relaxation spectrum 
and the pore-size distribution: First, the porous medium 
is pictured as an aggregation of isolated pores, which al- 
lows an unambiguous notion of the pore size if individual 
pores are of simple geometry. This may not necessar- 
ily require each pore to be closed. A periodic, symmet- 
ric arrangement of pore space [30] connected by narrow 
channels may well be considered as such, as long as the 
inter-pore diffusive flux either balances out or becomes 
negligible. Attempts to map between the pore size- and 
the surface-enhanced T2— distributions may work very 
well for systems such as mono-disperse bead packs, food 
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elements composed of suspended spherical voids such as 
cheese, and the class of sedimentary rocks such as clean 
sandstones. 

Second, in the absence of diffusive flux among such 
pores, the so-called fast- diffusion condition is met so 
that the relaxation spectrum is dominated by the slow- 
est mode for each pore the rate for which becomes 
proportional to the respective surface-to- volume ratio. 
For a simple isolated spherical pore of radius a, for 
example, the condition involves a single dimensionless 
parameter pi) 



Poa/D < 1. 



(1) 



The macroscopic parameter po characterizes the uniform 
depletion rate at the interface and is directly related to 
the probability of depletion p and enter the Robin's con- 
dition in the form of [Dh • V -t- po] ^(r) = on the bound- 
ary with h being the unit vector normal to the interface. 
When this condition prevails, it was noted [3T] that the 
surface-induced depletion rate Wg (or so-called surface- 
enhanced T2 relaxation rate [15] in the NMR context) is 
directly related to the surface (S)-to-volume(V) ratio of 
the pore: 



Po 



V 



(2) 



This simple relationship had been applied widely in ex- 
amples mentioned above^l^ US |M1 IHl [211 113 il [23 US 
l27l [28] . For classes of porous media with a broad vari- 
ation of its geometrical properties and strong diffusive 
coupling among its pore-constituents, these assumptions 
may break down. For example, in the oil-exploration, 
problems have been long recognized for the class of rocks 
in which the pore shape and the lithological composition 
of the matrix become complex. The complications in- 
duced by the heterogeneous, extended pore space in MR 
as well as other physical properties such as the electrical 
and hydraulic conductivity pose a fundamental challenge 



and invites active debates. As we will show in section III 



this condition is facilitated by the near-uniformity of the 
slowest eigenmode (See Eq 22 1 . When the pore geome- 
try has more than one length scale (Eqp^, the slowest 
mode acquires more pronounced spatial variation and the 
condition is relatively poorly met. 

The third assumption often made is that the surface- 
relaxation occurs with a uniform strength (i.e. p(r) — po) 
throughout the interface even though an inhomogeneous 
p is the norm, rather than an aberration, in natural me- 
dia. In porous rocks, there are several mechanisms for 
the surface-enhanced relaxation [311 [SB [35], and 

they often involve the stochastic distribution of magnetic 
minerals in the matrix or the local interfacial morphol- 
ogy. The strength of p arising from such origins would ac- 
quire inhomogeneity across the pore-grain interface, but 
microscopic details of such a variation is not known quan- 
titatively in general. 

Given the lack of such information, it is not entirely 
possible to dismiss the following observation: That for 




FIG. 1: An artistic rendition of situations considered in the 
text. Panel (a) shows a simple spherical pore with a uni- 
form Po (its strength indicated on the shell with a uniform 
gray color). Inside, the color represents the local population 
density (white:high, black:low) as evolved from a uniform ini- 
tial distribution in the long time limit. Panel (b) shows the 
same, but with a non-uniform p(r) on the shell. Panel (c) 
shows a more complex pore geometry but with uniform p. 
(d) shows the same with p(r) with a potential disruption on 
the registry between the local density and its pore-geometrical 
parameters. 



the population evolution of a collection of isolated pores 
of varying sizes, all satisfying the condition above (Eq[T]), 
it is possible to construct a collection of identical pores 
whose size amin is chosen such that PmaxO,min/D <^ 1, 
but with a distribution of p values G [pmim Pmax] as- 
signed to each, that will yield the identical evolution. It 
is worth noting that the question of whether there is a 
unique mapping between the eigenspectrum and a given 
boundary geometry has been posed in more abstract and 
stronger tcrmsjTl [SJ |3]. The hypothetical situation for 
an NMR relaxometry as posed above violates a weaker 
form of isospectral criterion as it involves the behavior of 
the overall population decay {A4{t)) obtained under both 
the uniform initial distribution and detection sensitivity 
profile (see Eq 21 in the following) under the condition 
K <C 1. Aside from mathematical rigor, questions arise 



as to whether the direct mapping between the geometry 
and the relaxation rates could remain useful for a general 
p(r) profile. Figure [l] summarizes the core issues in the 
form of the slowest mode profiles (casually rendered here 
for illustrative purpose) inside porous media with simple 
and moderately complex shapes (panels a and c). They 
are then further complicated by the presence of inhomo- 
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geneous p(r) which is incommensurate with the variation 
of the boundary shape (panels b and d). 

Several authors had considered the effect of an inho- 
mogeneous p(v). Wilkinson et al incorporated the in- 
homogeneous p in a toy model [S7] in reduced dimen- 
sions. Kansal and Torquato considered a numerical tech- 
nique to derive the effective trapping rate for a mix- 
ture of partially aborbing traps in the context of bio- 
logical systems [H]. Valfouskaya et aZ[38] considered a 
non-uniform absorption on randomly distributed sites in 
reconstructed porous media. Arns et al used numerical 
simulations with a focus on the cross-correlation between 
the relaxation spectrum and the transport property [39]. 
While the latter touches directly on one of the impor- 
tant practical issues, it provides little insight beyond the 
complications due to the pore-geometry issue (first and 
second conditions) alone. 



II. SETTING UP THE PROBLEM 

This work is concerned with consequences of allow- 
ing either of or all three assumptions above to break 
down. We aim to develop a method that incorporates 
the two components (geometry vs. inhomogeneous p) 
on an equal footing. To be precise, consider a pore 
space (V denotes its pore- volume in the following) de- 
fined by the solid(grain)-pore interface S, its area des- 
ignated as S) in a Euclidean space of dimension de- A 
physical property (such as polarized spin) is carried by 
molecules (or agents) diffusing through V with its mobil- 
ity characterized by the local diffusion tensor D{y). We 
allow such molecules to get absorbed (or killed) by a cer- 
tain mechanism at the boundary, if the property we are 
tracking is their population density, or allow the physi- 
cal property to be drained upon contacting the interface 
with a certain probability. The strength of such surface- 
localized depletion mechanism is controlled by a param- 



eter p(r), (See Eqs 10 and 27 1 which defines the coarse 



grained strength of absorption/depletion/relaxation. For 
an isotropic system, the probability of depletion per col- 
lision with the boundary, p, is related to p via|37l HOI HI] 
p{v) = f Y i-p/2 where the Brownian particle moves 
with continuous step sizes uniformly distributed in the in- 
terval [— e, e] for each of the directions during the time 
step. Inhomogeneity in p(v) may arise through spatial 
variation in p and/or D /e. The microscopic mechanism 
for the draining probability (p) will affect the texture, the 
spatial prifile, of p(r), but we will derive our main results 
without assuming a specific pattern for p(r). Using the 
standard bra-ket notation|33 US] the local population 
distribution at time t is represented in terms of the basis 
functions {|r >} as < rj^* >f (= ^'(r,t)), and overlap in- 
tegral between two such functions < '^\^' > is equivalent 
to dr^'*(r)^(r). The basis functions {|r >} satisfy the 
orthogonal property: < r'|r >= ^(r — r') where i5(r — r') 
is the (ig-dimensional Dirac-delta function with the nor- 
malization Jy i5(r — r')dr = 1, integrated over the pore 



volume V and r' G V. In the following section, we con- 
sider the diffusion equation according to which an initial 
state \^ >t=o evolves and consider the spectral property 
of the associated boundary value problem. Specifics of 
pore shape variation is incorporated into more generic 
spectral features of the modes. The varying degree of 
break down for the second and the third conditions is 
then systematically studied via the changes refiected on 
the spectra for a range of values in k (Eqs|l]and 15) and 
the dimensionless parameter 



< \Spir)\ > _< |p(r) - Pol > 



Po 



(3) 



with Po being the interfacial average of p(v). Obviously, 
< Sp{r) >— 0. What are the main observable conse- 
quences for allowing cr ^ in a natural porous media? 
In this work, we focus on the changes in the eigenvalue 
and the spatial mode profile of the slowest mode treat- 
ing a as the small parameter. We point out that this 
was largely motivated through mapping our Helmholtz 
problem to that of Schrodinger problem for the particle 
in a partially absorbing box in the imaginary time do- 
main and treating dp(r) as a perturbativc potential[44]. 
The mapping underscores the significance of the statis- 
tics and symmetries of the modes, especially the ground 
state, which are not as readily apparent in the tradi- 
tional Green's function formalism. This then prompts us 
to ask whether the effects should be more pronounced 
where the inhomogeneity is commensurate with the vari- 
ations in the underlying boundary shape, which in turn 
affects the spatial profile of the modes. As a corollary, 
it follows that the faster modes, having little spectral 
overlap with the spatial variation of p(r), unless p(r) is 
self-affine, may be less sensitive compared to their slower 
counterparts. These aspects should be considered on an 
equal footing along with the diffusive coupling H51 ITf] 
in affecting the spectral properties of the problem as we 
will elaborate in sections |III| and |IV[ For a class of exper- 
imental diffusion probes, parallels with the spectroscopy 
of a quantum particle [44] yield useful insight and has led 
to novel applications [m |49] even with a = 0. 

The organization of the rest of the paper is as follows: 
We consider in section |TTT] the case of a uniform depiction 
strength po on the boundary and offer an expanded ac- 
count of observations (Eqs 13|16 and 22 1 that had been 
made earlier [44]) on the spectrum for a general boundary 
geometry with uniform p. In section |IV[ we further de- 
velop for the spatially varying p(r), establishing a set of 
fundamental relationships linking the uniform and non- 
uniform cases through a perturbativc solution for the 
fractional changes of the eigenvalues and their weights. 
Eq[60] represents the central result of this perturbative 
approach. As a concrete example, we apply the method 
for a spherical pore with a general angular variation of 
p(r) in section [v] and obtain solutions for a specific binary 
(5p(r) texture. 



4 



III. UNIFORM p 

We start by considering the simpler case of a uniform 
p(t) = pq and a general local diffusion tensor D. Intro- 
ducing the flux operator J 

J = -D • V (4) 

and the Hamiltonian operator TL: 

7i = V • J = -V -D • V, (5) 

the evolution of ^' as dictated by continuity follows 

dtYl> >^ -n\'li > (6) 

which formally links the slowest depletion rate of our 
diffusion-relaxation problem with the ground state en- 
ergy of the analogous quantum mechanical system. [44] 
In real life MR relaxometry, diffusion of polarized spin- 
carrying molecules suffers an additional depletion in the 
bulk of the fluid if there exists a static field gradient. 
This dephasing may be eliminated via an experimental 
technique and therefore we neglect it for simplicity and 
consider only the depletion localized at the interface. 
Let us consider a partially absorbing boundary: 

< r|n(r) • Jl^f po < r|* > forr e S (7) 

where n(v) is the unit surface normal vector at the 
boundary point r pointing into the solid matrix. Time 
evolution of an initial distribution \^ >o can be ex- 
pressed as a linear superposition of the set of eigenmodes 
{|</>o>}(p = 0,l,2,...)of7^ 

OO 

Each eigenmode 10° > satisfies the equation 

n\^l >= A>o > (9) 
and the boundary condition at interface S: 

< r|n(r) • J|0° po < r|0° > forr e S. (10) 
Multiplying Eq|9]by < we obtain 

<^l\V-3\^l>^Xl (11) 

where the left hand side, upon inserting the complete 
set of basis functions / =|r> Jdr<r|, becomes the 
volume integration of 

V • (< 0°|r X r|J|0O >)- < r|J|<^o > -V < 0^|r > . 

(12) 

Combined with the boundary condition, we obtain the 
following expression for the eigenvalue: [44 

AO =« 0°|po|0° » + < 0°|J • ■ J\cf>l > (13) 



where << 



>>E 



. da. Employing the spatial 



representation, the right hand side is equivalent to 



po|0°(r)pdf7+ / i5„Mr)(Va0°(r))(V^0O(r))dr 



(14) 

where a, f3 = x,y,z. The result breaks the rate associ- 
ated with each eigenmode into two channels: a surface 
integral involving pQ and a volume integral involving spa- 
tial variaion of the mode, J|0p > and allows us to gen- 
eralize the criterion of slow- and fast-diffusion regime [5T| 
beyond the simple pore geometry. For the spherical pore, 
Brownstein and Tarr had shown that the dimensionless 
parameter k = pga/D controls the qualitatively distinct 
behavior for the time evolution of Based on the ob- 
servation that the distinction originates from the degree 
of spatial fluctuation in the slowest mode (See Eq 22 be- 
low) , we generalize the criterion by defining k parameter 
as the ratio of the two terms in Eq{T4] for the slowest 
mode, 100 >: 



•Jl 



> 



po /dr(V0g)^ 



« 0S]|po|0^ » 



(15) 



If the diffusive flux in the bulk dominates, (k — > oo), 
the slowest rate becomes independent of po, while in the 
opposite hmit, we have Aq po <f da{(j)'^)^ — poS/V. We 
therefore identify the length-scale parameter 



/rfr(V0°)2 



(16) 



for given po and D, as the relevant size that separates 
the distinct regimes for an arbitrary pore shape. It is 
interesting to note that £ is reminiscent of the A param- 
eter introduced by Johnson et al\50\ in the context of 
electrical conductivity in general porous media. 

EqfTS] may be also used to investigate the effect of 
changes in po and D as induced via control parameters 
such as the temperature, T. We obtain 

rfAO ^0|rfPom ^2 « 4>l\po\Scf>l »(17) 



dT 



dT 



< J • ■ J|0° > +2 < </)°|J • ■ J|J0O > 



where |(5</>p >= ^|0p >• Note that the slowest mode is 
the most sensitive to ^ in the so-called fast-diffusion 
limit, where the contribution from the surface-integral 
dominates over the second term in Eq|T4] 

Here we also prove that the spectral weight of the ex- 
cited modes in an initially uniform distribution is closely 
related to the spatial fluctuation of the slowest mode. To 
do so, let us take the initial state 



^-0(0 



^e(r.V) 



(18) 
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FIG. 2: A schematics of a typical evolution of a popula- 
tion M with a uniform boundary condition with po- At early 
times, the slope matches that of poS/V, and later times, it 
approaches that of the slowest eigenmode with Aq. (Each 
represented by broken curves). The spectral weight for the 
excited modes (W^) is also indicated. 

where Q{x) = 1 if the boolean condition x is satisfied, 
otherwise. Its subsequent time evolution follows 

l*>*=E^°e"'°l<> (19) 

with its spectral distribution given by 

5° >o= 4y / <(r)dr. (20) 

V V Jv 

The total population follows 

„ all 

M{t) = — / ^,{r)dr = ^ \s°fe-'^". (21) 

It can be shown that the fraction of the initial population 
belonging to the excited modes (i.e. g ^ 0) is 

W"^J2 \4\' = l-^oP >-\<^"o> n (22) 

which means that the total weight for all excited modes is 
directly proportional to the mean-square variance of the 
lowest mode. It is straightforward to show using Eqs. [6] 
and |7] that the slope at very early times, if the initial 
distribution is uniform, 

d_ ^ Jdr<r\n\^>t ^ § dapo < r|^ >t 
M{t)dt ^' j'dv<v\^>t /dr<r|*>f 

(23) 



which reduces to poy using the property that M.{t) is 
uniform in the limit i — > 0. Note that the result applies 
also to the case of disjoint pore systems if one under- 
stands S and V to represent the total interface area and 
the pore volume respectively. However, the range in time 
during which this assumption remains valid depends on 
the pore shape and is not a universal property as it de- 
pends on the boundary geometry of a given system. The 
following sum rule also follows: 

all ^ 

T.K?K = p4 (24) 
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and from the positive definiteness of A^'s, it also follows 
that Ag is bounded 

Po^ > A°. (25) 

The properties described in this section are schematically 
depicted in Figure [2j 



IV. NON-UNIFORM p 

In this section, we derive expressions for the changes in 
the eigenvalues and their modes when the boundary con- 
dition varies from point to point. For this, we should in- 
troduce another set of eigenmodes {4>p\ with eigenvalues 
{Ap} as opposed to the superscripted eigensystem {0°} 
and A^'s for uniform p(r) = pQ. Using the self-adjointed 
property of 7i, all eigenvalues are shown to be real, and 
their associated eigenmodes may be represented as real 
functions, as we choose to do so in the following. Each 
\(j)p > now satisfies 

n\(t)p >= Xp\(j)p > (26) 

and the non-uniform boundary condition 

< r|n(r) • J|(/)p >=< r|/9(r)|(/)p > forreE. (27) 

Following the steps that led to Eq{T3] we obtain 

Xp =« (j)p\p{r)\(Pp » + < (Pp\3 ■ D-^ ■ 3\(f>p > . (28) 

Our primary interest is now on the difference between the 
two eigensystems, without and with spatial variations in 
p, represented in terms of {|(/)° >}'s. Figure [s] shows the 
schematics of changed properties in comparison to the 
uniform po case, as we derive the details in this section. 

With the definition of 6p and a given in EqjS] We start 
by assuming that the eigenmode \4>p > of the inhomoge- 
neous case be expressed as a perturbation of the corre- 
sponding mode in the homogeneous counterpart, \<j)p >■ 

>= cp(|0° > +1% >). (29) 
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One may further decompose \S(j}p > into two components: 
|0P> = Cp(|0°>+^apg|</>°>) + [/-7']|.^p> 



all 

CpiY.apgK > +\Qp >] 



(30) 



where V = J^q I'/'? >< 'f'ql is a projection operator 
into the Hilbert space spanned by the set of eigenmodes 
{I't'q >}• O'pp = 1 by definition. The normahzation con- 
dition < (l)p\4>p >= 1 gives 



< QplQp >) 



-1/2 



(31) 



\Qp >, which by definition should satisfy < (f>q\Qp >— 
for all q's, represents the part oi \(j)p > that cannot be 
accounted for via a linear superposition of >'s. That 
\Qp > is not a null function, despite the fact that one 
can realize a least-square fit approximation [5T] with any 
given precision to an arbitrary function inside V, foll ows 
since the set {(pp} satisfies a boundary condition (Eq 10) 
while satisfies another (Eq27l. Strictly speakmg, 

any linear combination of 0p's, without the \Qp >, cannot 
satisfy the inhomogeneous boundary condition, while at 
the same time, the role of \Qp > is probably secondary 
to '^q^papq\4''l >, which needs to be verified. Therefore, 
in the following, we first focus on getting a^p's in terms 
of A°'s and the overlap integrals of i5jo(r) and 0^'s. We 
then evaluate \Qp > which satisfies 

—6p{r)\(j)p>-n-3\Qp>+pa\Qp>=0 (32) 

Cp 

on the boundary and is a solution to the inhomogeneous 
equation: 



(n-Xp)\Qp > 



E(^p 



A°)ap,|<^° > . (33) 



> 



To be systematic, we obtain an approximate \Qp 
with the perturbative solution Op^g's at the given stage 
((5p™ with m = 1,2, . . . ) substituting for the inhomoge- 
neous source term on the right hand side. It is assumed 



that IQp^m >^ IQp > as m 



oo. 



Substitute Eq|30] into Eq|26] and multiply both sides 

by <0°|, 



y^pUpq 



< 



b"q\n\Qp >-- 



XpQ,pq . 



(34) 



From the self-adjointedness of H and the boundary con- 
dition, we make the crucial observation: 



< 



(See Appendix-[A| where 

« «<>gl'5p|</>p »= <p 0°(r)(5p(r)(^p(r)da, 



(35) 



(36) 



0.1 



exp(- 1 X„ ) 



_L 



exp(- 1 X\) 
J 



0.0 0.2 0.4 0.6 0.8 

t 

FIG. 3: Schematics for the difference between the population 
evolution with uniform po and an inhomogeneous p(r) with 
the finial slopes given by Aq and Ao as indicated by the broken 
curves. W — represents the change in the spectral weight 
distribution. 



Using this, one can show 



(Ap - X")apq = (^ aprSpqr + 5pqp) 



(37) 



where we introduce the surface overlap integrals 6pqr and 



IP' 



^Pqr^ =« (t)l\5p\4>l » 



^P7py €\^P\Qp,^n » . 



(38) 
(39) 



It formally follows that 



«cl^"q\Sp\c^p»+Spq (40) 



which, using Eq 30 gives a recursive equation for a 



■pq- 



(1 - Sp,) 



s 



J-pq 



Xp — x^ 



I'l'^aprSpgr + Spqp^'^ + Spq. (41) 



Upon first iteration, and using the notation (5p™ associ- 
ated with the perturbative approximations \Qp,m >, we 
obtain 



^pq 



= ^pq + (1 - ^pq)[i 

Spqr Sprp 



''11 



Ap-AO Ap-AO^V 



(E 

r^p 



EE 

r^p s^p 



Xp — A^ Ap 

Spqr Sp, 



Spqr Spf.p 



E 



A„ — A^ Ap — A^ 



V)(y)^} 



7 



providing a way for a systematic expansion in powers of 
^ in a manner analogous to the diagrammatic expansion 
of a particle interacting with the perturbative potential. 
Truncating the iterations at second order, we obtain, 



where we used the orthogonality properties of the 
Qo and defined < (SA)^ >= ^ J^A^dr 



and < A >= ^ Adr where A stands for Qo, 



and 

„r ^0 



"q or syQ. 



E 

r^p 



Spqr 5 p. 



rp 



Xr, 



AO \ 



{^f + 0{5p^)] (42) 



since 6p^6pqr = 0{Qp^m) x Sp^ or higher. To obtain a 
formal solution in an algebraically closed form, we put 



p = q in Eqs 37 and obtain 

5 _ SppqSp, 
-PPy 2- AO -A, 



Ap — Ap — Spp 



E 



5 p. 



PI 



gp / ^ N 2 



^P\ 



S 



A,, 



V 



(43) 



where |(5(/)p >= J2q^p^pq\^q > +\Qp > and therefore the 
last term in Eq 43 contributes terms of order 0{dp^) and 



To make further progress beyond Eq|43] we need to 
solve for \Qp > and its surface integrals with respect to 
(5p(r)(/)0(r). Assuming that 10° > and A^ are known for 
po, obtaining \Qp >, apq and Ap in a mutually consistent 
manner constitutes the complete solution of the prob- 
lem which is not possible for a general pore geometry. 
Instead, we are interested in how changes in observable 
properties such as the slowest eigenvalue and its spec- 
tral distribution depend on the texture of 5p{r) and the 
spatial profile of the relevant eigenmodes that reflect the 
underlying boundary geometry via the k parameter in a 
perturbative scheme based on small a. 

As noted with Eq 42 and Eq|44] if we restrict our- 
selves to the second-order perturbation evaluation of i5Ap, 
we need only to construct |Qp,i > (i-e. \Qp > eval- 
uated using the first order perturbative solutions) as 
the contribution of Qp,m to 6Xp for m — 1,2,... is 



<< 



\Sp\Q: 



p,m 



©((5^™+^) or higher. We noted 
earlier that \Qp^i > arises from a distribution of source 
higher. Keeping only up to second order in 6p, and using that is the remnant of dp{r) as J2q '^pql'^'q > alone fails to 

^ account for its effect completely. Using Eqj33] and Eqj37] 

we obtain the following inhomogeneous Helmholtz equa- 
tion for (3p,i(r) to lowest order in a: 



the boundary condition for \Qp >, we transform Eq 
into the following alternative form: 



43 



\ \°-^n ^ V SPpqSPqP f'S 2 I.e. r \ 

Xp-X.^-dppp-^ (-) ^-[Spp-lpp) 

q^p q ^ 



where we define 



S 

(44) 



(H- Ap)|Qp,i >= |/p,i > 
where the source term is now given by 



(48) 



S. 



qp 



s 

V 



4P,{v){5p{v)fc[>l{v)da (45) 



l/p,i >= 



^ all 



I*; >« 0>/>i0; » ^ (49) 



and 



<j>^''^{v)5p{v)n-3Qp{v)da. 



(46) 



and take the p — q case. Note that Spp assumes the form 
of a self-interaction which, through the mixed-boundary 
condition, is mitigated by the Tpp term. We show in the 
following (section |v]) that, for the spherical pore geom- 
etry, these two terms largely cancel each other. It sug- 
gests that the reduction in Ap gained through the self- 
interaction term (Spp) is lost due to the rapid bending of 
the mode profile (Tpp) near the boundary. 

We also have the generalization of the result we ob- 
tained for the spectral weight for the excited modes 
(Eq|22|: 



By extending the surface localized function Sp{r) into a 
thin shell {dV^) of thickness e lining the interface in the 

pore space, and ensuring that e <C ^D/Xp, y^ZJ/A^, we 

define 



_ / eSp{r') forr € dV, (r' || (r - r'),r' e S) 

(50) 



q^O 



q^O 



V<dQl> -2V(< 0|] > ^ao, < > 



q^O 



< 00 >< Qo > + XI < <Pq >< Qo >) 



'^^^^'■^ { Otherwise 
to show that Eq|49]is equivalent to 

all 

- VSp\<l>l >, (51) 

the surface localized source function (5p(r)0O(r) projected 
onto the space spanned by the eigenmodes {|0p >}. In 
Appendix-[B| we show that \Qp > is given by the super- 
position of waves Gp emanating from the residual charge 



^ p,res • 



q^O 



\Qp >= (1 + GpSp) ^Gp\(Jp^res > 



(52) 



/(E"Og+<Qo|Qo>) 



(47) where Gp is the Green's function as defined in Eq B4 



For the second-order perturbation, the residual charge is 
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defined as 



\<Jp,res >= V5p\<j)l > -5pV-\<j)p > 



(53) 



Note that the last term is of order 0{Sp^) unless q = 
0. Therefore, taking only the q = contribution and 
rearranging, 



which is not necessarily limited to be on the interface. In 
the perturbative scheme, \(f)p > should now be replaced 
with J2q^pq\'l'q > +\Qp,i >■ The surface overlap inte- 



gral of \Qp > and Spqp (Eq 39 1 can now be put into the 
following form: 



E 



SPOqSpqO ,S 2 



E 



5 Pi 



Or 



AO - AO 



5proi^f + Oi6p').i59) 



Gp(ri,r2)(7'(r2,r3)<5p(r3)(/.^(r3) 



Mr3)^(r3,r2)-^(r2) (54) 

Cp I 

Closed solutions for Gp and the residual charge for a gen- 
eral pore shape and (5p(r) are not readily available, but 
evaluating Gp out of the basis set {(/)0}, and using or- 
thogonality of {0g}, we further obtain 



> 



AO - A 

q r^p q 



- « (l>°\Sp\(P° » apr (55) 
p 



which indicates that the requirement < (f>q\Qp >~ is 
obeyed in 0{Sp^). Qp{v) may be interpreted as the po- 
tential field induced by the residual charge distribution 
'^p,res and its contribution to the eigenvalue Ap is the in- 
teraction between the potential and the surface charge 
distribution associated with each mode, 5p{Y)(f)^q{Y). The 
potentail is subject to a destructive interference when 
<7p,res(r) has a rapid spatial fluctuation, and further 
weakened due to a veraging over the diffusion length 

in ^ 1 / w • (Eq B4 in the appendix) The leading order 



contribution of Qp to Ap, 5ppp of Eq |43[ is then 

X~ S ,01.1.0 «(tPq\5p\(tPp» 

Sppp^ ^ -l^X^apr « 4Pr\Sp\4>q» A • 

q r^p 1 P 



Finally, we arrive at a compact expression 



SXp = Sppp 



Sp 



9#P « 



— (-)' + 



5p{v)4Pp{v)Qp{v)da 
(57) 

for the change in eigenvalue valid up to second order in 
5p. For the slowest mode, we conjecture that |(3p=o >i 
as the surface source 5p4>^, is significantly weakened via 
convolution and commutation with V and the oscillatory 
kernel of wavelength ^Q. Taking p = 0, we obtain the 
fractional shift in the decay rate of the slowest mode: 



^0 

A8 



5poo S 
AO V 

5por 



SpOqSPqO 



EE^ 

q r^O 



AO- A 



oSprq 



SPqO _,S^^3 



A^V^ 



(58) 



Due to the presence of SXq in the denominator of the last 
term, \Qo > therefore makes a second order correction 
to the first term, only if Spoo ^ 0. Otherwise, its effect 
vanishes to second order in bp. We thus arrive at: 



(^Ao bpQQ S 



XI 



AS] V 



1 



SpOqSpqQ S 2 



SpOqSpqO 



EaO(AO-AO)^V 



-AO)^V 
OiSp') 



(60) 



The fractional change in the weight for the excited 
mode (Eq|47| can now be simplified: 



W 



EqMw^iv)' < > 



< > 



(61) 



noting that < Qq >~ 0, and < SQl >^< Qq >oc 
0{Sp'^). We also expect that the fiuctuations induced by 
Sp in the bulk of the pore when averaged over the pore 
volume will tend to vanish, < J2q O'0q4>q >~ 0; leaving 
the positive definite term above. (More precisely, one 
can show that these cancellations arise rigorously from 
the normalization condition, Eq 65 1 This shows that as 



the boundary condition becomes more inhomogeneous, 
the excited modes gains in weight in proportion to their 
mean-square fiuctuation, < {S(f>'^)^ >, but also weighted 
down by the l/(Ao— AO)^ factor and Sp^Q, overlap between 
(/)0 and 



0. The overall effect, however, is second order 
(56) in Sp at most (^ Ooq)- These results are schematically 
summarized in Figure [3j 

Following the steps taken for the uniform case, one can 
also show that the initial slope of population decay with 
a finite a satisfies: 



t^o Mit) dt ^ ' 



Pa 



S §Sp{r)^t{v)da 



V 



(62) 



As the second term vanishes for sufficiently short t when 
< rj^* >t is still uniform, it suggests that the initial slope 
may remain robust. Note however that the way the even- 
tual deviation off the initial slope sets in may be differ- 
ent from that of the uniform p^ case (even though the 
asymptotic value remains the same at poS/V) depending 
on how the depletion of population population proceeds. 
We also note that the sum-rule that we found for the 



uniform case (Eq 24 1 retains the same form, and can be 
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put into the form: 



E 



II 



SXg N 



1 + 2E 



r^q "-rq 



V<Ql> 



Pa 



S 
V 



(63) 



with SSq —< 



I* >o= E,.^ga<?rS°+ < Qgl^- >o. Here 
sjl is the overlap integral of \(f>'^ > with the initial distri- 
bution (defined earlier Eq 20 1 and \6(j)q > is as defined in 
Eq|30l Combining this witn Eq 24 and taking the lead- 
ing order in dp only, we have the following condition that 
has to be satisfied 



where jL{kL{n)r) is the spherical Bessel function asso- 
ciated with the angular momentum L. (Note that the 
modes defined thus are complex. However, the eigen- 
values are all real, and the results so far remain valid.) 
fci(n) denotes the radial functions each associated with 
the n— th nodal, positive solution of 



Po3L{aL{n)) + DkL{n)f^{aL{n)) = 



(67) 



where the dimensionless parameter a„ = kL{n)a and c„,l 
is the normalization constant so that /y 10° ^jyj(r)p(ir = 



1. Since / dnYLAjYj 



L'M' 



L.L'Sm,M' 



it follows [ 



(64) 



Note that it is in fact simply expressing the rigidity of 
the initial slope against the variational effect of 6p. To 
be complete, we note also that 



(65) 



which follows from the normalization requirement. 



C„,L 



V2aL{n) 



-3/2 



ai(n)2-(L+i)2 



jL{aL{n)) 



(68) 

In the analysis of Brownstein-TarrfSl' and its subsequent 
application to a variety of porous media, all modes with 
L 7^ are excluded from consideration because all rele- 
vant integrals vanish under the uniform boundary condi- 
tion and the isotropic initial state(L = 0,M = 0). Here, 
we consider the non-uniform 5p{v) parametrized via its 
overall strength a and the angular variation on the sphere 
of radius a, /(f^): 



SPHERICAL PORE 



5p{v) = (Tpaf{n) 



(69) 



As a solid example, let us consider the case of an iso- 
lated spherical pore. Although it is a three dimensional 
object, much of its MR response reduces to that of a 
one-dimensional system with a single controlling length 
scale. This is often overlooked and its properties have 
been casually interpreted as generic to pores with com- 
plex three dimensional morphology. To be more concrete, 
we sketch out analytic expressions and numerically eval- 
uate them for some of the properties for the spherical 
pore despite this reservation. As an intermediate step, 
extension of the methods developed here to non-spherical 
pore geometry |52j would be useful. Angular variations in 
the boundary condition bring in aspects of the extra di- 
mension and length scales more explicitly, (such as the 
terms contributing to Eq 30 with L > 0) and it would 



be interesting to study the effect on k (Eq{l5]) of (j}Q as it 
varies along the boundary of non-trivial geometry. Yet, 
to make an impact on SXq/Xq, their effects need to sur- 
vive the angular-averaging; furthermore, the associated 
eigenmodes should have sizable presence on the boundary 
of the pore. Therefore, one should approach with discre- 
tion the conclusions we draw from the spherical pores in 
the following. 

The eigenmodes for the uniform po inside the spheri- 
cal pore are separated into the radial part {jL{kr)) and 
the angular part (Ylm{(^ ,4'))- Instead of the generic in- 
dex g, we employ the set of indices {n,L,M) that char- 
acterize the eigenmode associated with the eigenvalue 
DkL{n)'^ in terms of its radial and angular parts: 



(r) = Cn.LjL{kL{n)r)YLM{^) 



(66) 



with which the modes with a finite angular momentum 
contribute to 5\q to further slow down the slowest mode. 
Note that, in a similar manner, these higher eigenmodes 
play increasingly significant role as the pore geometry 
further deviates and acquires more asymmetry and het- 
erogeneity. Due to the L = Q sy mm etry of the mode 

vanishes, and for 



>, the first order term in Eq 58 



the second order term, only states with the L, M— com- 
ponents that are present in the 5{v) profile contribute. 
Thus for the spherical pore, we have, up to second order 
in 5p, 



6\a _ SXa.a S\o,b 



-^0 



-^0 



-^0 



(70) 



where 



■5 An 



is the contribution from coupling to the 



eigenmodes {10" >}(q ^ p) and coming from IQo >• 
Let us examine the two second order contributions one 
by one. The first term, ^4if^, is 



-^0 



E/ ''l,0'-n. 



2 2 



Lk{ko{l)af jL{kL{n)af x 



(<^Po)^^o,L;M 
A^(A:i.L-AEj)4^- 



where we define 



(71) 



(72) 
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FIG. 4: An example of eigenvalues kL{n) for odd angular 
momenta L = 1, 3, . . . , 99 shown for n up to 200 for k = 
0.4161. The solid line (fci(l)a/7r ~ 0.51L°-^) is just a guide 
for the eyes. Shown in the inset is the angular factor Wo.iiM 
(Eq |72[ | for the hemi-spherical p(r) variation. 



Log JO [(5Ao,^,(n,L)/(5Ao,fl(l,l)] 
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FIG. 5: Contribution to 5Ao,a from each eigenmode with 
radial node index n = 1,2,... and angular momenum L{= 
1, 3, 5, . . .) for K = 0.416. The contour levels represent base-10 
logarithm of individual contribution S\o,a{n, L) normalized to 
the maximum value 5Ao,a(l, 1) ~ 0.0323. 




2000 4000 6000 8000 10000 

Number of included modes, N 




FIG. 6: Upper panel: Convergence of numerically evaluated 
(5Ao,a as total number of modes increases. All eigenmodes in 
the ranges 1 < n < 200 and 1 < L < 101 were found, and 
their individual contribution 5Ao,a(w, L) evaluated and sorted 
according to their magnitude. The graph shows SXo.aii) 
as the number of included modes N increases. Lower panels: 
Dependence of partially summed SX{n,L) on n{= 1,2,...) 
and L for k = 0.416. The left panel shows contribution from 
all modes with same n summed over L < 101. The right 
panel shows contribution with same L values summed over 
all n < 200. The solid lines are guides for the eyes with 
oc and oc L^^ respectively. 



as the Yl.m component in the harmonic expansion of 6p. 
Introducing the slowest rate for k ^ oo, 



(73) 



and using Eq|68| and aL{n) = ki^{n)a, we can put this 
into a form which displays its dependence on k explicitly 
for an arbitrary angular variation of 6p: 



S\o,a 



2 2 1 \ ^ 2 



0,L;M 



L,M,n 



Aq i^n.L ^ ^a) 



(k - \Y + alin) -{L+ i)2 - I) + alii) 
Let us consider a simple case with 



(74) 



(75) 



in which only L — odd modes with M = are present. 
Figure [4] shows the eigenvalues found for odd L up to 
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fco(l)n 



Cl,0 



jo(fco(l)a) 



SXg 1 

AS <t2 



0.01 
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1.4565 
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4.1614 
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41.6138 
208.07 



0.056176 
0.112001 
0.17599 
0.246325 
0.341266 
0.409545 
0.463496 
0.57835 
0.65411 
0.736232 
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0.976014 
0.995194 
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FIG. 7: Second-order contribution to the fractional shift of 
the slowest relaxation rate ((5Ao,a/Ao,o) for the hemispherical, 
binary distribution of p(r) for the spherical pore at various 
values of k. Up to 20000 modes were included to ensure good 
convergence for all k values. The contribution peaks around 
K ~ 2. 



10^, with n < 200 for each L. Figure [s] shows how much 
individual eigenmode with L, n contributes to SXq while 
Figure |6] shows the rapid convergence properties as the 
number of included mode increases. It also shows the 
partial contributions from all modes with a given L or n 
values. 

For small k ^ 1, Aq/Aoo oc k and Aoo/(A° ^ — Aq) 
becomes independent of p, and therefore one can see 
S^o.a/^o ~ while for k oo, we note that Aq/Aoo 
becomes independent of k, therefore 



A. 



^ K.L - AS! - 0.5)2 + ^^(„)2 _ + 0.5)2 ^ 

n 

so that SXq ^^/Xq oc k ^. Figure pj and Table fH show the 
numerically evaluated SXa^a/Xg for a wide range of k with 
the hemispherical Sp which bears out this observation. 
For a qualitative description, this may be roughly de- 
scribed as 



<JA^/A° 



2.3 



(77) 



The result suggests that, at least for the spherical pore, 
the system is most sensitive to the inhomogeneity in the 
intermediate range k £ [0.5, 10] in its second order con- 
tribution. Incidentally, this is where the two terms in 
the rate (as we noted earlier Eq 13 and also in [H]) 
are comparable and the system becomes most accom- 
modating of the perturbation 5p. If the former domi- 
nates (i.e. K ^ 1), it becomes too costly for the slow- 
est mode to deform itself from quasi-uniformity to ac- 
commodate Sp, while in the opposite case, the mode 



TABLE L Numeric values for ko, ci,o, jo{koa), 
for K = 0.01 to 200 obtained using summation up to 20000 
eigenmodes. 



amplitude near the boundary is severely reduced (i.e. 
jo{ko{l)a) , jL{kL(n)a) — > 0), and SXq becomes insensi- 
tive to a fractional change in Sp. 

More generally, for an eigenmode of angular variation 
with {L,M) to contribute at least a fraction /3 of SX/Xq, 
the strength of the corresponding component in the har- 
monic variation of Sp, loqx-.m would have to meet 



2 2 



LJo(fco(l)a)ji(/^L(?^)a 



AS) 



>f3 



which can be used to define the region of relevance in the 
{L, n}-plane for numerical evaluations. Even for pores 
without spherical symmetry, this criterion may be gen- 
eralized using the strength of the associated modes av- 
eraged over the interface that should replace and 
Ci oJo- It should be emphasized however, that, if the first 
order contribution survives, this second order effect may 
become overshadowed. 

Evaluation of the second term of Eq 70 ° 
involved for an arbitrary pore geometry, 
to evaluate first the V projection of the surface local- 
ized function Sp(r) and further its overlap integral with 
Gp^o{ri,r2), neither of each is available in a closed 
form. For spheres, however, we can clearly see from 
Eq |59| that it should make a vanishing contribution as 
<<^\Sp\(l)Q >>= as < r|0o >(x Yoo{i}). Comparisons 
to an exact solution [55 and numerical simulations [55] 
verify that 



A, 

as 



Ti— , IS more 
need 



we 



^Ao,j 



li j> do4>lSp{y)4>l{y) 







(79) 



for the second order contribution of Qo to SXq in spher- 
ical pores with an arbitrary Sp(v). At the same time, it 
is plausible that there exist boundary shapes for which 
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FIG. 8: Radial profile of the partially-projected, ^—averaged 
< 'PL,nm^P >e (the curves are multiplied by (~)^ 1 a 



for normalization) near the pore edge where L — 1 only and 
rim = kma/n as indicated for each curve up to 800 (see inset 
for blow-up for large rim)- The height of the normalized curves 
at r/a — 1 converges toward 0.75075 for k — 0.4 in this 
example. With inclusion of more modes with L > 1, the 
height approaches the value of 1.0 gradually as one would 
expect for a perfect representation. The con verg ence becomes 
progressively slower as k increases. (See Fig 10 1. 



(/>Q(r) develops a significant angular variation so that 
/ (/)q (5/900 dcr 7^ 0. In such a case, the first order contribu- 
tion in Eq |60| would dominate. \Qp > only contributes to 
its higher order modification. It is instructive to examine 
how the projection V and the residual source |tTo,res > 
behave in more detail. 

Let us first consider the numerical evaluation of the 



projected Sp, Vf{^) of the binary distribution (Eq 75 1 



For a delta-profile for Sp(r) in the radial direction with 
g^(5p(r) = at the boundary, VSp, as evaluated nu- 
merically with a large number of (/>q's, may approach 
Sp with an arbitrary precision jSlj. and yet fail to meet 
the zero-slope condition since all {(/>q}'s have a finite 
slope on the boundary unless po = 0. This discrep- 
ancy may hardly impact the accuracy of the surface in- 
tegrals / da(l)gSp{r)cj)Q in practice. Discrepancy between 
Sp and its projection VSp may become pronounced along 
the (de ~ 2) dimensional manifold (i.e. the equatorial 
line 9 — tt/2 in the hemispherical example) across which 
sharp changes in Sp occur at length scales smaller than 
£p. How this translates into an enhanced contribution 
to SXo^b can only be addressed numerically for a general 
Sp{r) texture. In the following, we investigate how the 
numeric VSp(r) representation behaves as a series-sum 
over a finite number of modes for the simple spherical 
model. 

Figure |8] shows how the radial delta-function like pro- 
file is approached with progressively larger number of ra- 
dial modes (with cutoffs kmax as indicated) averaged over 




FIG. 9: Strength of the projected 'PL,nm^P around 6 = ii/2 
at r = a across which it should go from —1.0 to 1.0 in a 
stepwise fashion. The curves show progressive refinement as 
we increase the number of included modes by increasing the 
maximum L up to 399 ( kL{nm)a/n was fixed at 1000 for each 
included L). 



6 G [0, 7r/2] . The oscillatiotory tail is due to the finite cut- 
off . The inset shows details near the boundary for larger 
cutoff values. With only L = 1 modes included, the value 
on the boundary converges to the value of ~ 0.75, signif- 
icantly short of 1.0. Inclusion of L > 1 modes remedies 
this, but its convergence is significantly impeded as k 
increases. 

Figure|9]shows how the angular profile (Eqj75 1 is repro- 
duced with progressively larger number of angular modes 
(with cut-off L,nax as indicated). The radial cutoff is set 
at kmaxtt/TT — 1000. In the hemispherical case, a moder- 
ate value of Lmax seems sufficient to achieve an accept- 
able convergence, although we observe that its rate slows 
down as k increases. 

We monitored the following dimensionless parameter 
as a measure of convergence for the 7^-projected Sp to 
the actual Sp (as an overlap integral with (j)p{r) over the 
pore- volume) in contributing to SXpj,: 

§da2Q{Sp{v^))Sp{v^)ct^l{v^) 

(80) 

where Q{x) is the heavy side step function, — 1 for a; > 
and ~ otherwise. Its convergence is largely determined 
by the spectral weight of modes {</'q^p(r)} present in 
(5(r)(/)p(r). To be systematic, we define the partial pro- 
jection operator VL,n^ = I]"" \4'n,L,o >< 0n,L,o| that 
projects onto a subspace spanned by the first rim ra- 
dial modes for each L. Panels (a) and (b) of Figure 
\T0\ show the angular average of the partial projection 
X^l""" 'PL,n„,Sp\(j)p > for values of Lmax and (equiv- 
alent to krnax)- In panel (b), small but rapid oscillations 
observed in radial direction are immaterial as they are 
averaged out when convoluted over the pore volume. As 
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the convergence of Cp is largely controlled by the spec- 
tral compostion of the source profile, f{^), this may no 
longer hold for a complex Sp texture. [57] 

For the hemi-spherical Sp, Cp for p ~ becomes 



0=0 = lim 



fe<Ti/2 



d(T2 § rfcTs < r2|^L,n„|r3 > 



Te<TT/2 



(81) 

where is the radial width of the projected delta peak 
on the boundary for the given n™, and is given by the 
quarter of the wavelength associated with t he m ode with 
khinm)- Em = | fc^(n„) - Panel (c) of Figure 10 shows this 
Cp=o as we progressively increase the number of modes in 
its numerical evaluation for both large and small values 
of K. It is worth noting that while the convergence of Cp 
(and therefore ^ Ao,b) is quite slow, the series sum for SXq ^ 
is much more rapid due to the factor of .p ^ .q (in Fig- 

urc p| . In a numerical simulation that employs random 
walkers with a fixed step size, one would be effectively 
truncating the series summation at a wavelength compa- 
rable to the stepsize. The apparent strength of a in such 
simulations may then deviate from what corresponds to 
a fully converged Cp in the figure. Our result provides a 
guide on how one may correct for such artifacts in a sys- 
tematic manner. Numerical simulations employing large 
number of random walkers, which uses continuous step 
sizes to alleviate such an issue, are underway for various 
types of pores. 
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FIG. 10: Convergence of the projected 5p(j)^ to the actual. 
Panel (a) shows that the value of the partial projection (L = 1 
modes only) with k = 200, averaged over the hemispherical 
shell approaches ~ 0.746, falling short of the expected 1.0 
even after summing over 40000 radial eigenmodes. Panel (b), 
which shows the radial profile of the projected (5p<^o for = 
0.4, indicates that good convergence is achieved when modes 
with L up to 15 or more are included. Panel (c) shows how the 
convergence Cp=o 1.0 slows down as k increases. With k = 
200, even after including up to 3 x 10® modes, the convergence 
is not quite complete. 



We considered the consequence of the spatially varying 
boundary condition [D{Y)n ■ V -|- p(r)]^'(r, t) = for the 
spatio-temporal evolution of the local density \E'(r,t) of 
an attribute carried by diffusing entities. It has direct 
relevance on the local magnetic (polarization) density of 
fluid molecules in the magnetic resonance relaxometry 
widely used for various porous media for their charac- 
terization. We examined the spectral properties of the 
governing Helmholtz equation and their relationship to 
the boundary geometry and the texture of its control- 
ling parameter p(r) . Using only the general properties of 
the modes and the boundary conditions they satisfy, we 
showed that each eigenvalue can be expressed as a sum 
of two parts (EqpSj) : one is in the form of a surface in- 
tegral directly involving p{v), the other being a volume 
integral which involves the diffusive flux of the mode. 
The direct relationship between the slowest eigenvalue 
and the surface-to-volume ratio of the pore is recovered 
when the flrst term dominates over the second, and we 
derived the generalized parameter k (Eq |15[ ) that quanti- 
fies the boundary between distinct regimes with observ- 
able consequences in the evolution of an initial distribu- 
tion. We also showed that the weight of the slowest decay 
mode in the overall relaxation of the attribute is dimin- 
ished in direct proportion to the rms spatial fluctuation 



of the mode (Eq 22 1 . Traditionally, the weight for all 



modes other than the slowest is often interpreted as rep- 
resenting small pores, a notion increasingly invalid as the 
pores become extended through diffusive coupling and 
acquire complex geometry. We clarified issues regard- 
ing the time-domain evolution of ^'(r, t) which originate 
from such inadequate interpretation of the spectral distri- 
bution, (Eq 20 1. Building on this, we then introduced 
spatially varying p(y) and obtained the perturbative so- 
lution in a =< \Sp\ > /pq. The results show how the 
effect of 5p{r) manifests itself as the shift in the lowest 
eigenvalue and is controlled by overlap integrals of Sp{r) 
with the associated mode. (Eq 36) We also showed and 
verified numerically that the initial slope of the overall de- 
pletion remains robust (See Eq64l. These results were 



derived without making any specific assumption about 
the pore geometry, relying only on the self-adjointedness 
and boundary conditions of the problem. We show that 
the first order contribution vanishes when the base sys- 
tem has symmetry so that the overlap integral with the 
mode (j)p, § (/)p(r)(5p(r)0p(r)dcr — 0. When the boundary 
geometry varies in a complex manner, the slowest mode 
00 (r) itself acquires significant spatial variaion even un- 
der the uniform po, and the incommensuracy between 
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Sp{r) and 0q's may result in a significant non-zero first or- 
der contribution (Eq |57| , which overshadows the second- 
order effect. In the opposite hmit, where the texture of 
Sp{r) is such that its variation occurs on length scales 
much shorter than the pore geometrical length-scale £ 
(Eq |16[ ) , the effect of such finely inhomogeneous Sp should 
be muted via diffusive averaging-out. This is the case 
considered in Valfouskaya et al |38j which considered a 
texture with a random variation uncorrelated beyond the 
voxel size, much smaller than the typical grain size, in a 
stochastically generated 3D porous medium. Our numer- 
ical simulations on random glass bead packs with similar 
textures of Sp also yielded results consistent with these 
observations. |41] Towards the opposite limit, we applied 
out theory to a case where the impact of a finite a may be 
most pronounced (section^. Extensive numerical anal- 
ysis including up to 3 x 10^^ eigenmodes was performed 
for the simple case of a spherical pore for a wide range of 
K = poa/D. We propose that the fractional change in the 
slowest eigenvalue is the most effective probe into varia- 
tions in the boundary condition. We examined how much 
each of the eigenmodes contributes to the change in the 
slowest rate (Figure [s]) , and obtained the overall second 
order contribution for a wide range of K (Fig [7]) .The lat- 
ter is observed to peak around the value of k ~ 2.0, and 
follows roughly cx k and cx 1/k at either end. Our result 
provides a useful theoretical framework and quantitative 
bounds for more complex situations addressed mainly 
through numerical simulations. [Ml ESI |4r Further de- 
velopment through comparison to exact solution [55 and 
systematic numerical simulations |56j are underway. 
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APPENDIX A 



Here we show that 

< (^qiniS^bp >-< <5</)p|H|0° >= c;^ « ^°\6p\cj)p » 

(Al) 

where > is an eigenmode with the uniform boundary 
condition, \(f)p > is an eigenmode with the inhomogeneous 
boundary condition, Cp, its normalization constant with 
>= Cp(|0° > +\Sct)p >) and 7^ = V • J with J = -D • 
V. For notational brevity, we choose the representation 
in which all the eigenmodes are real functions. We start 
by putting the first term on left hand side as 

^ (V • (/)i;(r)J%(r) - (J(50p(r)) • V(/)i;(r))dr (A2) 

and the second term as 

/ (v.<5(/.p(r)J(^^(r)-(J</.0(r)).V%(r))dr (A3) 



and using the Gauss's theorem to put the left hand side 
into: 



^°{r)h ■ 3S^p{r) - 5<f>p{r)n ■ J(/.°(r)) da. (A4) 



Substituting S(j)p{r) = Cp-^0p(r) — 0p(r), the integrand 
becomes 

</>°(r)(n . - n . J^°(r)) _ ( _ ^O(r))^ . j^0(,) 



which upon using the boundary conditions turns Eq |A4| 
into: 



j^4>,{^)5p{v)4>p{v)da. (A5) 



APPENDIX B 



Here we derive the particular and homogeneous solutions 
for \Qp > at a given stage in the perturbative iteration 
and show how Qp is related to 5p and |(/>p >. We start by 
noting that Qp{y) is the solution to the inhomogeneous 
Hclmholtz equation 



{n-\p)\Qp >-|/p> 



(Bl) 



for the \Qp > function where Ap is the perturbative eigen- 
value consistent with the boundary condition as satisfied 
by the solution from the previous stage, /p(r) is given ac- 
cordingly in Eq |33| The boundary condition that should 
be satisfied by Qp is 

-Sp{r) J2 ap«'/'g(r) - "(r) • JQp(r) + p(r)Qp(r) = 0. 

^ Q 

(B2) 

Using the projection operator V — J2q l*^!? 'I'ql O'^ito 
the Hilbert space spanned by >, \Qp > is related to 



> via 



\Qp >=-{l-V)\^p> 



(B3) 



from which it follows that VlQp >— 0. Our aim is to seek 
a formal solution for Qp(r) using the Green's function 
approach. Consider the following Green's function 



(7^-Ap)Gp(r,ri) = (5(r-ri). 



(B4) 



We consider its representation in the basis functions of 
an eigen-system {\^q >} that satisfies 

{n-eg)\^q>^0 (B5) 

and the general condition 

(po - n{r) ■ J)e,(r) = C(r)e,(r) (B6) 
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on the boundary with the function C(r) to be chosen for 
convenience. Choice of C(r) = —Sp{i^) amounts to solv- 
ing for {(f>q}, while the lowest order perturbation would 
amount to the choice of C(r) = 0. 

The Green's function Gp may now be given in terms 
of >: 



G,(r,r')=<r|^|^, > 



AS 



< e,|r' > (B7) 



and it satisfies the boundary condition 

(po - n{r) ■ 3)Gp{r, n) = C(r)Gp(r, ri). (B8) 

To obtain \Qp > as a perturbative solution for the source 
function fp and Gp, we multiply Eq Bl by Gp(ri,r) from 
the left side and integrate over r to get 

GpHlQp > -XpGplQp Gplfp > . (B9) 

Using Stoke's theorem, this becomes 

inGp)\Qp> + GpUJ\Qp>^i3G)UQp>-XpGp\Qp> 
= Gplfp > (BIO) 

where |e indicates a surface integral. Replacing {TiGp) — 
Xp Gp + / and using the boundary conditions for Qp 



(EqB2l and Gp (EqB8l, after rearranging, we get 

(/ - GpUSp{r) + C(r))) \Qp Gp\ap,res > (BU) 
where we introduce the residual source: 



The general solution for \Qp > admits possible addi- 
tion of the solution, \Qp >, for the homogeneous coun- 
terpart of Eq 33 We now show that \Qp >— on physi- 
cal grounds. First, note that the homogeneous equation 
for \Qp > and the original problem for \(j)p > becomes 
identical. The same applies to the boundary condition. 
This suggests we may take the homogeneous solution 
\Qp > to be identical to \(j)p > up to a constant factor a, 
\Qp >— a|(/)p >, with a to be determined by consistency 
requirements: 



y^p > > MQp >) = \Qp 



and the boundary condition 



> 



(B13) 



[-n • J + po]{\Ql > +\Q^ >) + Spq)p^O (B14) 



Bll 



Since 



where |Qp > is the particular solution of Eq " 

> term should vanish by its construction, and re- 
placing \Qp > with a\(j)p > on the left hand side turns 
the boundary condition into 



(1 -I- a)Sp{r) 0p(r) = Jp(r) ^p{r) (B15) 



which therefore requires a = 0, i.e. we should take 

\Q^ 0. 
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